# split plot design
sp <- read.table('ecol 563/splityield.txt', header=T)
sp[1:8,]
 
# reorder the levels of density from low to high
levels(sp$density)
sp$density <- factor(sp$density, levels=c('low','medium','high'))
levels(sp$density)
 
# aov uses the classical approach to fitting split plot designs
# requires specifying an error term that shows the hierarchical structure
out0 <- aov(yield~irrigation*density*fertilizer + Error(block/irrigation/density/fertilizer), data=sp)
summary(out0)
#the lowest level is not needed
out0 <- aov(yield~irrigation*density*fertilizer + Error(block/irrigation/density), data=sp)
 
# examine the two-factor interactions
interaction.plot(sp$fertilizer, sp$density, sp$yield)
interaction.plot(sp$fertilizer ,sp$irrigation, sp$yield)
interaction.plot(sp$density, sp$irrigation, sp$yield)
 
# split plot design as a mixed effects model
 
# use nlme package
library(nlme)
out0a <- lme(yield~irrigation*density*fertilizer, random=~1|block/irrigation/density, data=sp)
anova(out0a)
 
# use lme4 package
detach(package:nlme)
library(lme4)
out0b <- lmer(yield~irrigation*density*fertilizer + (1|block/irrigation/density), data=sp)
anova(out0b)
 
# alternatively we can create a variable the identifies the split plots and split split plots
 
# split plot
sp$sp1 <- paste(sp$block, sp$irrigation, sep='.')
sp[1:18,]
# split split plot
sp$sp2 <- paste(sp$block, sp$irrigation,sp$density, sep='.')
sp[1:18,]
 
# now each plot is uniquely identified. To use lmer requires new notation
# does not work
out0c <- lmer(yield~irrigation*density*fertilizer+(1|block/sp1/sp2), data=sp)
# this works
out0c <- lmer(yield~irrigation*density*fertilizer+(1|block)+(1|sp1)+(1|sp2), data=sp)
anova(out0c)
 
detach(package:lme4)
library(nlme)
 
# with nlme we can use either variables in the same notation
out0a <- lme(yield~irrigation*density*fertilizer, random=~1|block/sp1/sp2, data=sp)
# drop nonsignificant interactions
out1a <- update(out0a, .~.-irrigation:density:fertilizer-density:fertilizer)
anova(out1a)
 
# display final model
library(lattice)
# show structure of the design including the blocks
xyplot(yield~density|fertilizer+irrigation, data=sp, groups=block, type='o')
# same graph as previous graph but written as a pannel function
xyplot(yield~density|fertilizer+irrigation, data=sp, groups=block, panel=function(x,y,groups,...){
panel.superpose(x, y, groups, type='o', ...)})
 
# obtain predicted mean (marginal) from model and add to data frame
sp$pred <- predict(out1a, level=0)
sp[1:8,]
 
# first attempt
xyplot(yield~density|fertilizer+irrigation, data=sp, groups=block, panel=function(x,y,subscripts,groups,...){
# change color of observed profiles
panel.superpose(x, y, subscripts, groups, type='o', col='grey70', ...)
# this does not work--we are plotting too much
panel.xyplot(x,sp$pred[subscripts], col=2, pch=8, type='o')})
sp[1:20,]
 
# redo graph using only first 3 means from data frame in each panel
xyplot(yield~density|fertilizer+irrigation, data=sp, groups=block, panel=function(x,y,subscripts,groups,...){
panel.superpose(x, y, subscripts, groups, type='o', col='grey70', ...)
panel.xyplot(x[1:3], sp$pred[subscripts][1:3], col=2, pch=8, type='o')})
 
# save lattice graph as object and use latticeExtra pkg to place strips on top and left
library(latticeExtra)
xyplot(yield~density|fertilizer+irrigation, data=sp, groups=block, panel=function(x,y,subscripts,groups,...){
panel.superpose(x, y, subscripts, groups, type='o', col='grey70', ...)
panel.xyplot(x[1:3], sp$pred[subscripts][1:3], col=2, pch=8, type='o')}) -> my.graph
useOuterStrips(my.graph)
 
# compare graph to output of summary table: match estimates to what is shown in graph
summary(out1a)$tTable
round(summary(out1a)$tTable,4)
 
# repeat summary graph interchanging fertilizer and density to examine second significant interaction
 
xyplot(yield~fertilizer|density+irrigation, data=sp, groups=block, panel=function(x,y,subscripts,groups,...){
panel.superpose(x, y, subscripts, groups, type='o', col='grey70',...)
panel.xyplot(x[1:3], sp$pred[subscripts][1:3], col=2, pch=8, type='o')})
 
# the density values are connected in the wrong order--define factor order to match data frame 
sp[1:8,]
sp$fertilizer <- factor(sp$fertilizer, levels=c('N','P','NP'))
levels(sp$fertilizer)
 
# redo graph and save it as an object
xyplot(yield~fertilizer|density+irrigation, data=sp, groups=block, panel=function(x,y,subscripts,groups,...){
panel.superpose(x, y, subscripts, groups, type='o', col='grey70',...)
panel.xyplot(x[1:3], sp$pred[subscripts][1:3], col=2, pch=8, type='o')}) -> my.graph2
useOuterStrips(my.graph2)
 
#### repeated measures example ######
 
turtles <- read.table('ecol 563/turtles.txt', header=T)
turtles[1:8,]
dim(turtles)
 
# plot data showing repeated measure structure
xyplot(plasma~treat|sex, groups=subject, type='o', data=turtles)
 
# reorder treatments so levels are in time order
turtles$treat <- factor(turtles$treat, levels=c('fed','fast10','fast20'))
xyplot(plasma~treat|sex, groups=subject, type='o', data=turtles)
 
# fit full interaction model accounting for repeated measures structure
turtle1 <- lme(plasma~treat*sex, random=~1|subject, data=turtles)
anova(turtle1)
 
# drop nonsignificant interaction
turtle2 <- lme(plasma~treat+sex, random=~1|subject, data=turtles)
anova(turtle2)
 
# the sex effect is not significant--this does not match the earlier graph
# refit model using lmer and carry out mcmc
detach(package:nlme)
library(lme4)
turtle2a <- lmer(plasma~treat+sex+(1|subject), data=turtles)
out.mc <- mcmcsamp(turtle2a, n=10000)
 
# 95% credible interval for sex does not include zero
# Bayesian approach disagrees with frequentist approach
HPDinterval(out.mc)
fixef(turtle2a)